function elements = calculateLimitStage(elements,model)
ne = model.countElements();
for i=1:ne
    e = model.getElement(i-1);
    f = e.getFace();
    if elements(i).isCrack==0 && elements(i).isActivated==1
        %% check if the crack was initialized
        [center,normal,isCrack] = computeCrack(e,elements(i));
        if isCrack == 1
            elements(i).isCrack = 1;
            disp(['warning: crack initiated at element ' num2str(i)]);
            
            %% realize positive domain
            [PosDomain,NegDomain] = checkPositiveDomain(f,center,normal);
            elements(i).PosDomain = PosDomain;
            elements(i).normal = normal;
            elements(i).tangent = [-normal(2) normal(1)]';
            
            %% compute crack length
            [crackLength,X1,X2] = ...
                computeCrackLength(f,center,normal,PosDomain,NegDomain);
            elements(i).crackLength = crackLength;
            elements(i).X1 = X1;
            elements(i).X2 = X2;
            
            %% compute element's area
            elements(i).area = f.computeVolume();
            
            %% compute Kdxi & Kxid
            elements(i).Kdxi = zeros(8,2);
            elements(i).Kxid = zeros(2,8);
            elements(i).Kxixi_1 = zeros(2,2);
            ng = e.getNumGaussPoint();
            B = zeros(3,8);
            t = elements(i).thickness;
            
            for j=1:ng
                xi = e.getGaussPoint(j-1);
                w = e.getGaussWeight(j-1);
                jac = f.jacXAt(xi);
                detJac = det(jac);
                
                B = e.computeB(f,B,xi);
                C = elements(i).constitutiveMatrix;
                
                Gc = computeGc(elements(i),B);
                elements(i).Kdxi = elements(i).Kdxi + B'*C*Gc*(w*detJac*t);
                
                Ge = computeGe(elements(i),xi);
                elements(i).Kxid = elements(i).Kxid + Ge'*C*B*(w*detJac*t);
                
                elements(i).Kxixi_1 = elements(i).Kxixi_1 + ...
                    Ge'*C*Gc*(w*detJac*t);
            end
            
            %% compute Kxixi
            J = computeJ(elements(i));
            Ctilde = computeCtilde(elements(i));
            elements(i).Kxixi = elements(i).Kxixi_1 + ...
                crackLength*J'*Ctilde*J*t;

            %% initialize renh and fully open state
            elements(i).renh = [0 0]';
            elements(i).active = 1;
            elements(i).active_old = 1;
            
%             elements(i).Kdxi
%             elements(i).Kxid
%             elements(i).Kxixi
        end
    end
end
